Antje Hofmann: Learning Syntactic Relations with the Mole Task

RePsychLing in SMLP2022

Author

Reinhold Kliegl

Published

August 24, 2022

1 Background

Children’s (age: 4-8 years)reaction times in a task teaching them syntactic relations.

1.1 Overview

  • Original analysis is by Antje Hofmann.
  • MixedModels.jl version
  • Addition of new chunks illustrate
    • selection of parsimonious LMM using random-effects PCA
    • plotting conditional means
    • illustration of borrowing strength

2 Readme

2.1 Design

2.2 Variables

  • Subj: Participant ID (renamed from ID; random factor)
  • Item: Word ID (random factor)
  • age: 4 - 8 years
  • Block (within-Subj/within-Item):
    • 1st Learning
    • 2nd Learning
    • Disruption
    • Recovery
  • Target(renamend fom targetness)
    • non-syllable target
    • syllable target
  • rt: response time

3 Setup

3.1 Environment

3.2 Packages

First attach the MixedModels.jl package and other packages for plotting. The CairoMakie.jl package allows the Makie graphics system [@Danisch2021] to generate high quality static images. Activate that package with the SVG (Scalable Vector Graphics) backend.

Code
using AlgebraOfGraphics
using Arrow
using CairoMakie       # graphics back-end
using CategoricalArrays
using Chain
using DataFrames
using DataFrameMacros  # simplified dplyr-like data wrangling 
using MixedModels
using MixedModelsMakie # diagnostic plots
using ProgressMeter
using Random           # random number generators
using RCall            # call R from Julia
using StatsModels

using AlgebraOfGraphics: boxplot
using AlgebraOfGraphics: density

using MixedModelsMakie: qqnorm
using MixedModelsMakie: ridgeplot
using MixedModelsMakie: scatter
using MixedModelsMakie: caterpillar

ProgressMeter.ijulia_behavior(:clear);
CairoMakie.activate!(; type="svg");
  • The data are available as an arrow file.
  • Most of preprocessing was done with R in RStudio (see Hofmann_Maulwurf.Rmd).
  • Order of factor levels should be checked.
Code
dat = DataFrame(Arrow.Table("./data/Hofmann_Maulwurf_rt.arrow"))
transform!(dat, 
    :Target => categorical => :Target,
    :Block  => categorical => :Block,
    :age => (x -> x .- 6) => :a1, # center age (linear) at 6 years
    :rt => (x -> log.(x)) => :lrt)
describe(dat)

8 rows × 7 columns

variablemeanminmedianmaxnmissingeltype
SymbolUnion…AnyUnion…AnyInt64Type
1Subjp01p850Union{Missing, String}
2Item01PA.ogg32SG.ogg0Union{Missing, String}
3age6.210813.76.318.040Union{Missing, Float64}
4Block1st LearningRecovery0Union{Missing, CategoricalValue{String, UInt32}}
5TargetNon-target syllableTarget syllable0Union{Missing, CategoricalValue{String, UInt32}}
6rt2335.39228.02100.08389.00Union{Missing, Float64}
7a10.210813-2.30.312.040Float64
8lrt7.673675.429357.649699.034680Float64
  • Centering age at six years yields an interpretable GM
  • Factor levels can also be set when contrasts are defined (see below).
  • BoxCox check showed that reaction time rt [ms] should be transformed to speed [1/s] = [Hz]
  • Indicator variables for Target and Block generated in R.

4 LMM analysis

4.1 Contrasts

Code
contrasts = merge(
      Dict(:Target => EffectsCoding()),
      Dict(:Block => SeqDiffCoding()),
      Dict(nm => Grouping() for nm in (:Subj, :Item))
   );

4.2 Varying only intercepts LMM m_voi

Code
f_voi1 = @formula(lrt ~  1 + Block * Target * a1 + (1 | Subj) + (1 | Item));
m_voi1 = fit(MixedModel, f_voi1, dat; contrasts)
Minimizing 53    Time: 0:00:00 (15.47 ms/it)
Est. SE z p σ_Item σ_Subj
(Intercept) 7.6895 0.0371 207.05 <1e-99 0.0239 0.2637
Block: 2nd Learning -0.0920 0.0104 -8.87 <1e-18
Block: Disruption 0.0516 0.0126 4.09 <1e-04
Block: Recovery -0.0477 0.0125 -3.81 0.0001
Target: Target syllable -0.0108 0.0040 -2.69 0.0072
a1 -0.0896 0.0324 -2.77 0.0056
Block: 2nd Learning & Target: Target syllable -0.0057 0.0104 -0.55 0.5854
Block: Disruption & Target: Target syllable -0.0262 0.0124 -2.11 0.0347
Block: Recovery & Target: Target syllable 0.0020 0.0123 0.16 0.8741
Block: 2nd Learning & a1 0.0347 0.0091 3.83 0.0001
Block: Disruption & a1 0.0214 0.0108 1.97 0.0486
Block: Recovery & a1 -0.0159 0.0108 -1.47 0.1409
Target: Target syllable & a1 -0.0052 0.0035 -1.46 0.1430
Block: 2nd Learning & Target: Target syllable & a1 -0.0113 0.0090 -1.25 0.2111
Block: Disruption & Target: Target syllable & a1 0.0113 0.0108 1.05 0.2933
Block: Recovery & Target: Target syllable & a1 0.0086 0.0107 0.80 0.4248
Residual 0.3004

4.3 Extract indicator variables

Code
X = modelmatrix(m_voi1)
dat.trng = X[:,2];
dat.drpt = X[:,3];
dat.rcvr = X[:,4];
dat.trgt = X[:,5];
describe(dat)

12 rows × 7 columns

variablemeanminmedianmaxnmissingeltype
SymbolUnion…AnyUnion…AnyInt64Type
1Subjp01p850Union{Missing, String}
2Item01PA.ogg32SG.ogg0Union{Missing, String}
3age6.210813.76.318.040Union{Missing, Float64}
4Block1st LearningRecovery0Union{Missing, CategoricalValue{String, UInt32}}
5TargetNon-target syllableTarget syllable0Union{Missing, CategoricalValue{String, UInt32}}
6rt2335.39228.02100.08389.00Union{Missing, Float64}
7a10.210813-2.30.312.040Float64
8lrt7.673675.429357.649699.034680Float64
9trng-0.014562-0.750.250.250Float64
10drpt-0.0528019-0.5-0.50.50Float64
11rcvr0.0524467-0.25-0.250.750Float64
12trgt0.00331492-1.01.01.00Float64

Switch to indicator variables and refit the model.

Code
f_voi2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1 + (1 | Subj) + (1 | Item));
m_voi2 = fit(MixedModel, f_voi2, dat; contrasts)
Est. SE z p σ_Item σ_Subj
(Intercept) 7.6895 0.0371 207.05 <1e-99 0.0239 0.2637
trng -0.0920 0.0104 -8.87 <1e-18
drpt 0.0516 0.0126 4.09 <1e-04
rcvr -0.0477 0.0125 -3.81 0.0001
trgt -0.0108 0.0040 -2.69 0.0072
a1 -0.0896 0.0324 -2.77 0.0056
trng & trgt -0.0057 0.0104 -0.55 0.5854
drpt & trgt -0.0262 0.0124 -2.11 0.0347
rcvr & trgt 0.0020 0.0123 0.16 0.8741
trng & a1 0.0347 0.0091 3.83 0.0001
drpt & a1 0.0214 0.0108 1.97 0.0486
rcvr & a1 -0.0159 0.0108 -1.47 0.1409
trgt & a1 -0.0052 0.0035 -1.46 0.1430
trng & trgt & a1 -0.0113 0.0090 -1.25 0.2111
drpt & trgt & a1 0.0113 0.0108 1.05 0.2933
rcvr & trgt & a1 0.0086 0.0107 0.80 0.4248
Residual 0.3004

They are equivalent.

4.4 A zero-correlation parameter LMM m_zcp

Code
f_zcp1 = @formula(lrt ~  1 +  Block * Target * a1 + zerocorr(1 + Block * Target | Subj) + (1 + a1 | Item));
m_zcp1 = fit(MixedModel, f_zcp1, dat; contrasts);

show(issingular(m_zcp1))
VarCorr(m_zcp1)
Minimizing 365   Time: 0:00:00 ( 2.40 ms/it)
  objective:  2393.633425634436
true
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0742531 0.2724942
Block: 2nd Learning 0.0362822 0.1904790 .
Block: Disruption 0.0156401 0.1250605 . .
Block: Recovery 0.0133038 0.1153419 . . .
Target: Target syllable 0.0002450 0.0156539 . . . .
Block: 2nd Learning & Target: Target syllable 0.0000000 0.0000000 . . . . .
Block: Disruption & Target: Target syllable 0.0032559 0.0570605 . . . . . .
Block: Recovery & Target: Target syllable 0.0006594 0.0256785 . . . . . . .
Item (Intercept) 0.0009502 0.0308250
a1 0.0006752 0.0259853 -0.83
Residual 0.0767542 0.2770456

Again, check the equivalence.

Code
f_zcp2 = @formula(lrt ~  1 + (trng+drpt+rcvr) * trgt * a1  + 
                zerocorr(1 + (trng+drpt+rcvr) * trgt | Subj) + (1 + a1 | Item));
m_zcp2 = fit(MixedModel, f_zcp2, dat; contrasts);

show(issingular(m_zcp2))
VarCorr(m_zcp2)
Minimizing 365   Time: 0:00:00 ( 1.51 ms/it)
  objective:  2393.633425634436
true
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0742531 0.2724942
trng 0.0362822 0.1904790 .
drpt 0.0156401 0.1250605 . .
rcvr 0.0133038 0.1153419 . . .
trgt 0.0002450 0.0156539 . . . .
trng & trgt 0.0000000 0.0000000 . . . . .
drpt & trgt 0.0032559 0.0570605 . . . . . .
rcvr & trgt 0.0006594 0.0256785 . . . . . . .
Item (Intercept) 0.0009502 0.0308250
a1 0.0006752 0.0259853 -0.83
Residual 0.0767542 0.2770456

4.5 A complex parameter LMM m_cpx

Code
m_cpx = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + trgt * (trng+drpt+rcvr) | Subj) + (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_cpx))    # not ok
show(m_cpx.PCA[:Subj])     # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_cpx)) 
VarCorr(m_cpx)
Minimizing 1328      Time: 0:00:02 ( 1.63 ms/it)
true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .      .
 trgt         -1.0    1.0     .      .      .      .      .      .
 trng          0.35  -0.35   1.0     .      .      .      .      .
 drpt          0.06  -0.06   0.35   1.0     .      .      .      .
 rcvr         -0.42   0.42  -0.23  -0.48   1.0     .      .      .
 trgt & trng   0.26  -0.26  -0.12   0.09  -0.16   1.0     .      .
 trgt & drpt  -0.4    0.4   -0.07   0.14  -0.01  -0.76   1.0     .
 trgt & rcvr   0.61  -0.61   0.16   0.04  -0.3    0.48  -0.85   1.0

Normalized cumulative variances:
[0.441, 0.667, 0.8072, 0.906, 0.9581, 1.0, 1.0, 1.0]

Component loadings
 
               PC1    PC2    PC3    PC4    PC5    PC6    PC7    PC8
 (Intercept)  -0.46   0.14   0.36  -0.2   -0.28   0.12   0.05  -0.71
 trgt          0.46  -0.14  -0.36   0.2    0.28  -0.12  -0.05  -0.71
 trng         -0.19   0.43   0.17   0.76  -0.0   -0.4   -0.13  -0.0
 drpt         -0.1    0.49  -0.59   0.15  -0.19   0.57   0.13   0.0
 rcvr          0.26  -0.42   0.31   0.47  -0.39   0.51  -0.15  -0.0
 trgt & trng  -0.31  -0.36  -0.49  -0.01  -0.52  -0.34  -0.38  -0.0
 trgt & drpt   0.39   0.44   0.14  -0.29  -0.16   0.02  -0.72  -0.0
 trgt & rcvr  -0.46  -0.18  -0.04   0.08   0.6    0.33  -0.53  -0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
──────────────────────────────────────────────────
     model-dof  -2 logLik       χ²  χ²-dof  P(>χ²)
──────────────────────────────────────────────────
[1]         28  2393.6334                         
[2]         56  2353.2482  40.3853      28  0.0611
──────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0757935 0.2753062
trgt 0.0001016 0.0100809 -1.00
trng 0.0334971 0.1830221 +0.35 -0.35
drpt 0.0198200 0.1407834 +0.06 -0.06 +0.35
rcvr 0.0191718 0.1384624 -0.42 +0.42 -0.23 -0.48
trgt & trng 0.0020619 0.0454079 +0.26 -0.26 -0.12 +0.09 -0.16
trgt & drpt 0.0145508 0.1206268 -0.40 +0.40 -0.07 +0.14 -0.01 -0.76
trgt & rcvr 0.0085501 0.0924670 +0.61 -0.61 +0.16 +0.04 -0.30 +0.48 -0.85
Item (Intercept) 0.0000000 0.0000000
a1 0.0004752 0.0217982 +NaN
Residual 0.0768486 0.2772158

The deviance improves, but we end up with an overparameterized LMM. Remove the two small VCs of the three Block x Target interactions (i.e., the third contrast of the Block factor).

4.6 A parsimonious parameter LMM m_prm

We remove one of the VC for trgt * trng contrast interaction, that is one of the three interaction terms.

Code
m_prm1 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + trgt *(drpt + rcvr) + trng  | Subj) + 
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm1))  # not ok
show(m_prm1.PCA[:Subj])   # not ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm1, m_cpx)) 
VarCorr(m_prm1)
Minimizing 1329      Time: 0:00:02 ( 1.84 ms/it)
true
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .      .
 trgt         -1.0    1.0     .      .      .      .      .
 drpt          0.04  -0.04   1.0     .      .      .      .
 rcvr         -0.41   0.41  -0.47   1.0     .      .      .
 trgt & drpt  -0.43   0.43   0.22  -0.04   1.0     .      .
 trgt & rcvr   0.63  -0.63   0.02  -0.31  -0.88   1.0     .
 trng          0.34  -0.34   0.33  -0.23  -0.12   0.18   1.0

Normalized cumulative variances:
[0.4716, 0.7125, 0.8294, 0.9405, 0.9938, 1.0, 1.0]

Component loadings
                PC1    PC2    PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.5    0.01  -0.27   0.34   0.24   0.07  -0.71
 trgt          0.5   -0.01   0.27  -0.34  -0.24  -0.07  -0.71
 drpt         -0.08   0.63   0.22  -0.41   0.61   0.09  -0.0
 rcvr          0.28  -0.48   0.39   0.32   0.64  -0.18   0.0
 trgt & drpt   0.36   0.45  -0.35   0.35   0.04  -0.65   0.0
 trgt & rcvr
  -0.47  -0.23   0.12  -0.42  -0.04  -0.73   0.0
 trng         -0.24   0.34   0.72   0.44  -0.33  -0.07   0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + drpt + rcvr + trgt & drpt + trgt & rcvr + trng | Subj) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.6334                          
[2]         48  2338.3644   55.2691      20  <1e-04
[3]         56  2353.2482  -14.8838       8     NaN
───────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0759040 0.2755068
trgt 0.0000989 0.0099437 -1.00
drpt 0.0193595 0.1391386 +0.04 -0.04
rcvr 0.0188340 0.1372370 -0.41 +0.41 -0.47
trgt & drpt 0.0106015 0.1029636 -0.43 +0.43 +0.22 -0.04
trgt & rcvr 0.0081492 0.0902726 +0.63 -0.63 +0.02 -0.31 -0.88
trng 0.0349529 0.1869571 +0.34 -0.34 +0.33 -0.23 -0.12 +0.18
Item (Intercept) 0.0008589 0.0293075
a1 0.0006867 0.0262043 -0.91
Residual 0.0763357 0.2762891

We don’t lose goodness of fit, but are still overparameterized. We remove CPs for Target.

Code
m_prm2 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) + 
                  zerocorr(0 + trgt  | Subj) +
                          (1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm2))  #  ok
show(m_prm2.PCA[:Subj])   #  ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm2,  m_cpx)) 
VarCorr(m_prm2)
Minimizing 1001      Time: 0:00:02 ( 2.01 ms/it)
false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.22   1.0     .      .     .
 drpt & trgt  -0.36   0.26  -0.1   -0.07   1.0     .     .
 rcvr & trgt   0.59  -0.01  -0.27   0.12  -0.87   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.3441, 0.5974, 0.7402, 0.8568, 0.9478, 0.9939, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.07  0.0   -0.04   0.61  -0.58   0.15
 drpt         -0.09   0.61  0.0   -0.05  -0.62  -0.47   0.11
 rcvr          0.3   -0.49  0.0    0.52  -0.16  -0.58  -0.17
 trng         -0.26   0.36  0.0    0.83   0.09   0.3   -0.08
 drpt & trgt   0.47   0.44  0.0   -0.1    0.37  -0.14  -0.65
 rcvr & trgt  -0.59  -0.22  0.0   -0.14  -0.27   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.6334                          
[2]         42  2341.3448   52.2886      14  <1e-05
[3]         56  2353.2482  -11.9034      14     NaN
───────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0757927 0.2753047
drpt 0.0193698 0.1391753 +0.04
rcvr 0.0188467 0.1372834 -0.41 -0.47
trng 0.0351278 0.1874242 +0.33 +0.34 -0.22
drpt & trgt 0.0099421 0.0997100 -0.36 +0.26 -0.10 -0.07
rcvr & trgt 0.0066351 0.0814561 +0.59 -0.01 -0.27 +0.12 -0.87
trgt 0.0002651 0.0162812 . . . . . .
Item (Intercept) 0.0008262 0.0287433
a1 0.0006740 0.0259607 -0.93
Residual 0.0762062 0.2760546

Check the Item-related CP. It is very large, might be spurious.

Code
m_prm3 = let
    form = @formula(lrt ~  1 + trgt * (trng+drpt+rcvr) * a1  + 
                          (1 + drpt + rcvr + trng + drpt&trgt + rcvr&trgt  | Subj) + 
                  zerocorr(0 + trgt  | Subj) +
                  zerocorr(1 + a1 | Item));
    fit(MixedModel, form, dat; contrasts);
  end;

show(issingular(m_prm3))  # ok
show(m_prm3.PCA[:Subj])   # ok
show(MixedModels.likelihoodratiotest(m_zcp2, m_prm3, m_prm2, m_cpx)) 
VarCorr(m_prm3)
Minimizing 1129      Time: 0:00:02 ( 2.00 ms/it)
false
Principal components based on correlation matrix
 (Intercept)   1.0     .      .      .      .      .     .
 drpt          0.04   1.0     .      .      .      .     .
 rcvr         -0.41  -0.47   1.0     .      .      .     .
 trng          0.33   0.34  -0.21   1.0     .      .     .
 drpt & trgt  -0.37   0.25  -0.08  -0.07   1.0     .     .
 rcvr & trgt   0.61   0.04  -0.31   0.13  -0.86   1.0    .
 trgt          0.0    0.0    0.0    0.0    0.0    0.0   1.0

Normalized cumulative variances:
[0.35, 0.5979, 0.7408, 0.8585, 0.9488, 0.9943, 1.0]

Component loadings
                PC1    PC2   PC3    PC4    PC5    PC6    PC7
 (Intercept)  -0.52   0.03  0.0   -0.01   0.61  -0.57   0.16
 drpt         -0.12   0.62  0.0   -0.06  -0.61  -0.47   0.12
 rcvr          0.33  -0.47  0.0    0.52  -0.19  -0.58  -0.18
 trng         -0.27   0.36  0.0    0.83   0.08   0.31  -0.08
 drpt & trgt   0.45   0.47  0.0   -0.1    0.38  -0.15  -0.63
 rcvr & trgt  -0.58  -0.22  0.0   -0.15  -0.26   0.03  -0.72
 trgt          0.0    0.0   1.0    0.0    0.0    0.0    0.0
Model Formulae
1: lrt ~ 1 + trng + drpt + rcvr + trgt + a1 + trng & trgt + drpt & trgt + rcvr & trgt + trng & a1 + drpt & a1 + rcvr & a1 + trgt & a1 + trng & trgt & a1 + drpt & trgt & a1 + rcvr & trgt & a1 + MixedModels.ZeroCorr((1 + trng + drpt + rcvr + trgt + trng & trgt + drpt & trgt + rcvr & trgt | Subj)) + (1 + a1 | Item)
2: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + MixedModels.ZeroCorr((1 + a1 | Item))
3: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + drpt + rcvr + trng + drpt & trgt + rcvr & trgt | Subj) + MixedModels.ZeroCorr((0 + trgt | Subj)) + (1 + a1 | Item)
4: lrt ~ 1 + trgt + trng + drpt + rcvr + a1 + trgt & trng + trgt & drpt + trgt & rcvr + trgt & a1 + trng & a1 + drpt & a1 + rcvr & a1 + trgt & trng & a1 + trgt & drpt & a1 + trgt & rcvr & a1 + (1 + trgt + trng + drpt + rcvr + trgt & trng + trgt & drpt + trgt & rcvr | Subj) + (1 + a1 | Item)
───────────────────────────────────────────────────
     model-dof  -2 logLik        χ²  χ²-dof  P(>χ²)
───────────────────────────────────────────────────
[1]         28  2393.6334                          
[2]         41  2356.2006   37.4328      13  0.0004
[3]         42  2341.3448   14.8558       1  0.0001
[4]         56  2353.2482  -11.9034      14     NaN
───────────────────────────────────────────────────
Column Variance Std.Dev Corr.
Subj (Intercept) 0.0760499 0.2757714
drpt 0.0193875 0.1392389 +0.04
rcvr 0.0189413 0.1376274 -0.41 -0.47
trng 0.0350733 0.1872787 +0.33 +0.34 -0.21
drpt & trgt 0.0091677 0.0957481 -0.37 +0.25 -0.08 -0.07
rcvr & trgt 0.0059753 0.0773002 +0.61 +0.04 -0.31 +0.13 -0.86
trgt 0.0002813 0.0167715 . . . . . .
Item (Intercept) 0.0006681 0.0258480
a1 0.0005356 0.0231424 .
Residual 0.0763232 0.2762665

Perhaps not. We stay with LMM m_prm2

4.7 Compare models with goodness-of-fit statistics.

Code
let mods = [m_zcp2, m_prm2, m_prm1, m_cpx];
 DataFrame(;
    model=[:m_zcp2, :m_prm2, :m_prm1, :m_cpx],
    pars=dof.(mods),
    geomdof=round.(Int, (sum  leverage).(mods)),
    AIC=round.(Int, aic.(mods)),
    AICc=round.(Int, aicc.(mods)),
    BIC=round.(Int, bic.(mods)),
  )
end

4 rows × 6 columns

modelparsgeomdofAICAICcBIC
SymbolInt64Int64Int64Int64Int64
1m_zcp228290245024502639
2m_prm242287242524262709
3m_prm148276243424352759
4m_cpx56259246524662843

LMM m_prm2 is a defensible solution according to \(\Delta\) AIC, \(\Delta\) BIC suggests we should not bother with CPs.

Code
coeftable(m_prm2)
Coef. Std. Error z Pr(>
(Intercept) 7.68255 0.0387548 198.23 <1e-99
trgt -0.00948222 0.00442195 -2.14 0.0320
trng -0.0841465 0.0279289 -3.01 0.0026
drpt 0.0428818 0.0228927 1.87 0.0610
rcvr -0.0370297 0.0226366 -1.64 0.1019
a1 -0.0922379 0.0338516 -2.72 0.0064
trgt & trng -0.00176262 0.00963315 -0.18 0.8548
trgt & drpt -0.019929 0.01816 -1.10 0.2725
trgt & rcvr -0.00411596 0.0162158 -0.25 0.7996
trgt & a1 -0.00874271 0.0038894 -2.25 0.0246
trng & a1 0.0414053 0.0244065 1.70 0.0898
drpt & a1 0.0287421 0.0199004 1.44 0.1487
rcvr & a1 -0.0208157 0.0197058 -1.06 0.2908
trgt & trng & a1 -0.00234843 0.00840933 -0.28 0.7800
trgt & drpt & a1 0.0171553 0.0158677 1.08 0.2796
trgt & rcvr & a1 0.00435249 0.0141867 0.31 0.7590

4.8 Pruning fixed-effects

There is an option of pruning some higher-order interactions.

4.9 Diagnostic plots

Various visualizations are used to check whether or not data are defensibly modeled with an LMM. They may lead to removal of outliers, transformations of the dependent variable, and deliver valuable heuristic information to be followed up with exploratory post-hoc analyses or ideally replication of new insights gained this way. In practice, it appears that only severe violations will stop people from reporting a model.

4.9.1 Residuals over fitted

Code
scatter(fitted(m_prm2), residuals(m_prm2))

Looks like we missed some fast response times.

4.9.2 Q-Q plot

Code
qqnorm(m_prm2; qqline=:none)

Hm?

4.9.3 Residual distributions: observed vs. theoretical

Curves for residulas based on observed and theoretical values should correspond.

Code
let
  n = nrow(dat)
  dat_rz = (;
    value=vcat(residuals(m_prm2) ./ std(residuals(m_prm2)), randn(n)),
    curve=repeat(["residual", "normal"]; inner=n),
  )
  draw(
    data(dat_rz) *
    mapping(:value; color=:curve) *
    density(; bandwidth=0.1);
  )
end

Figure 1: Kernel density plot of the standardized residuals for model m1 versus a standard normal

They are a bit too narrow.

4.10 Conditional means of random effects

4.10.2 Borrowing-strength plots

Shrinkage refers to the adjustment of subject-level or item-level predictions by taking population estimates into account. The further a subject’s/item’s estimate is from the fixed effect or the more variable or less reliable the subject’s/item’s estimate, the more the prediction will be shrunk towards the population estimate. Alternative terms for shrinkage are “borrowing strength” (Tukey) and regularization. My favorite is actually Tukey’s because indeed we borrow strength from the population estimates to make predictions for individual subjects’ effects. The goal of this section to illustrate the results of borrowing strength.

Subject-related conditional means of random effects revealed information about individual differences beyond fixed effects. Would these results also be visible in unconditional means, that is when we compute GM and experimental effects within subjects (i.e., as fixed effects) without borrowing strength from the population estimates?

In the following plots, effect estimates based on alone on each subject’s data (i.e., no pooling of data, no borrowing of strength) are plotted in pink and the subjects’ conditional means shown in the caterpillar plots are plotted in blue. The arrows indicate how much a subject’s prediction is changed by borrowing strength from knowledge of the population estimates.

Code
shrinkageplot!(Figure(; resolution=(1000, 1200)), m_prm2)

Figure 3: Shrinkage plots of the subject random effects in model m_prm2

There is an outlier subject

Code
cm = raneftables(m_prm2);   
sort(DataFrame(cm.Subj), ["(Intercept)", :trng, "drpt & trgt"])
sort(DataFrame(cm.Subj), :trng)

52 rows × 8 columns (omitted printing of 1 columns)

Subj(Intercept)drptrcvrtrngdrpt & trgtrcvr & trgt
StringFloat64Float64Float64Float64Float64Float64
1p10-0.665206-0.08002240.176762-0.5029420.325148-0.281849
2p14-0.5171140.14949-0.104271-0.3804460.0865683-0.0397183
3p760.0792766-0.1005910.104864-0.309562-0.02850420.0205178
4p02-0.621993-0.1969990.260501-0.2516720.0704512-0.126543
5p360.421384-0.03653870.0872904-0.231847-0.1690420.145319
6p390.291843-0.1308840.143621-0.200882-0.01923620.00243074
7p74-0.32949-0.06072740.091318-0.1793690.0671451-0.0858905
8p670.03832130.00115633-0.0194392-0.1550160.0331431-0.0117934
9p06-0.2376990.1393920.026263-0.117284-0.06679530.0402144
10p540.313189-0.1620190.0466067-0.116629-0.03677660.0411771
11p23-0.0288494-0.07110220.0526879-0.103847-0.1270080.0803724
12p42-0.09738960.288231-0.166534-0.103720.121435-0.0338403
13p590.4216630.0211065-0.255664-0.10291-0.01234780.0649439
14p45-0.363685-0.2957540.169737-0.0976150.0435265-0.0933103
15p30-0.329959-0.149979-0.0897719-0.0879097-0.0134295-0.00936962
16p83-0.230431-0.0531665-0.160467-0.07026520.0132008-0.00841336
17p84-0.389883-0.2028250.444848-0.0491713-0.1993720.0205808
18p720.176311-0.0685062-0.0523226-0.0387344-0.1194970.118113
19p12-0.148148-0.04483020.0523053-0.0284382-0.06717190.0361489
20p040.138247-0.04483740.0391417-0.0278217-0.003231170.00442783
21p770.134277-0.09654760.0522975-0.0217883-0.08010990.0500009
22p79-0.0769189-0.01929190.0300986-0.0144482-0.03216940.010373
23p690.104371-0.0392821-0.0829739-0.01442770.02267180.0351996
24p55-0.198404-0.0636382-0.0785629-0.00935968-0.06243770.024083
25p270.0253214-0.00349694-0.00408308-0.00611823-0.09766170.0736573
26p280.150439-0.00601184-0.1515260.007032250.0683167-0.0248102
27p250.2676250.137029-0.04744690.00970517-0.03567270.0451016
28p640.2683220.002286270.03393510.0132941-0.01882420.0245811
29p570.0028259-0.00255108-0.01532360.01560190.0551938-0.0361024
30p53-0.07836780.02572590.02034980.01802720.0430869-0.0320897

It appears to be Subj “p10”.

5 Appendix

Code
versioninfo()
Julia Version 1.8.0
Commit 5544a0fab76 (2022-08-17 13:38 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin21.4.0)
  CPU: 12 × Intel(R) Core(TM) i9-8950HK CPU @ 2.90GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, skylake)
  Threads: 1 on 6 virtual cores